﻿#include "imgspec.h"
#include <Eigen/Dense>
#include <opencv2/core/eigen.hpp>

// Problem 1: Downsample and blur
cv::Mat ImageProcessor::downSampleImage(const std::string& imagePath,const int kernel) {
    cv::Mat img = cv::imread(imagePath, cv::IMREAD_COLOR); int k = kernel;
    if (k % 2 == 0)k++;
    cv::GaussianBlur(img, img, cv::Size(k, k), k);
    cv::resize(img, img, cv::Size(100, 100));
    return img;
}

// 构建相似图
Eigen::SparseMatrix<double> ImageProcessor::Laplacian_mat(const cv::Mat& inimg, const int rate) {
    int count = 0;
    cv::Mat img = inimg.clone();
    const int W = inimg.rows;
    const int H = inimg.cols;
    img.convertTo(img, CV_64FC3);
    img /= 255.0;
    std::vector<Eigen::Triplet<double> > triplets;
    Eigen::VectorXd D(W*H);

    const int diagC = 0;
    double posrate = rate * 1.0 / 10.0;
    int kernel = 3;
    cv::Vec3d sum = img.at<cv::Vec3d>(0, 0); int cnt = 0;
    for (int k = W - kernel - 1; k <= W + kernel; k++) {
        for (int l = H - kernel - 1; l <= H + kernel; l++) {
            int kk = k % W, ll = l % H;
            sum += img.at<cv::Vec3d>(kk, ll); cnt++;
        }
    }
    sum /= cnt + 1;
    for (int i = 0; i < W; i++) {
        for (int j = 0; j < H; j++) {
            double d = 0;
            if(!KNN)
            for (int k = W - kernel-1;k <= W + kernel; k++) {
                for (int l = H - kernel-1;l <= H + kernel; l++) {
                    int kk = k%W, ll = l%H;
                    cv::Vec3d cor=img.at<cv::Vec3d>(i, j);
                    if (diagC)cor = (img.at<cv::Vec3d>(i, j) + sum) / 2;
                    cv::Vec3d diff = cor - img.at<cv::Vec3d>(kk, ll);
                    double dist = diff.dot(diff) + pow((i - kk)*posrate / W, 2) + pow((j - ll)*posrate / H, 2);
                    dist = exp(-4 * dist);
                    //if (dist < 0.04)continue;
                    triplets.emplace_back(W * i + j, H * kk + ll, dist);
                    ++count; d += dist;
                }
            }
            if(KNN)
            for (int k = i - kernel; k <= i + kernel; k++) {
                for (int l = j - kernel; l <= j + kernel; l++) {
                    int kk = (k+W) % W, ll = (l+H) % H;
                    cv::Vec3d diff = img.at<cv::Vec3d>(i, j) - img.at<cv::Vec3d>(kk, ll);
                    double dist = diff.dot(diff) + (pow((i - k)*posrate / W, 2) + pow((j - l)*posrate / H, 2))*20.0;
                    dist = exp(-4 * dist);
                    if (dist < 0.04)continue;
                    triplets.emplace_back(W * i + j, H * kk + ll, dist);
                    ++count; d += dist;
                }
            }
            D[W * i + j] = d > 0?pow(d, -0.5):0;
           
        }
    }

    std::cout << "Total edges:" << count << '\n';
    std::cout << "average degrees:" << ((double)count / (double)(W * H)) << '\n';
    Eigen::SparseMatrix<double> A(W * H, W * H);
    A.setFromTriplets(triplets.begin(), triplets.end());A.makeCompressed();
    Eigen::SparseMatrix<double> Diag(D.asDiagonal());
    Eigen::SparseMatrix<double> I(W * H, W * H);I.setIdentity();
    Eigen::SparseMatrix<double> L = I + Diag * A * Diag;
    L.makeCompressed();
    

    return L;
}


// 幂迭代法求最大特征向量
Eigen::VectorXd ImageProcessor::computeSpecVector(Eigen::SparseMatrix<double>&L,int k) {
    const int n = 10000;
    //k *= 20;
    Eigen::VectorXd x = Eigen::VectorXd::Random(L.cols()) + Eigen::VectorXd::Ones(L.cols());
    x /= 2;
    for (int iter = 0; iter < k; ++iter) {
        x = L * x;
        if(iter %100==0)x.normalize();
    }
    x.normalize();
    return x;
}

// k-means聚类
void ImageProcessor::spectralClustering(Eigen::VectorXd& vec, const cv::Mat& img,
    cv::Mat& output) {
    int n = img.rows * img.cols; vec *= n;
    cv::Mat data(n, 1, CV_32F);
    for (int i = 0; i < n; ++i) data.at<float>(i) = vec(i);

    cv::Mat labels, centers;
    cv::kmeans(data, 2, labels, cv::TermCriteria(cv::TermCriteria::MAX_ITER, 200, 0.1),
        3, cv::KMEANS_PP_CENTERS, centers);

    output = cv::Mat::zeros(img.size(), CV_8UC3);
    for (int i = 0; i < n; ++i) {
        cv::Vec3b color = labels.at<int>(i) ? cv::Vec3b(0xcc, 0x77, 0x66) : cv::Vec3b(0xff, 0xff, 0xff);
        output.at<cv::Vec3b>(i / img.cols, i % img.cols) = color;
    }
}

// 可视化特征向量（色温映射：蓝色→绿色）
cv::Mat ImageProcessor::visualizeSpecVector(const Eigen::VectorXd& Spec, const cv::Mat& img) {
    int n = img.rows * img.cols;
    cv::Mat color_map(img.size(), CV_8UC3);

    // 归一化特征值到 [0, 1]
    double min_val, max_val;
    cv::Mat eigenMat;
    cv::eigen2cv(Spec, eigenMat);
    cv::minMaxLoc(eigenMat, &min_val, &max_val);
    double norm_val = max_val - min_val;
    if (norm_val < 1e-8) norm_val = 1.0; // 避免除零

    for (int i = 0; i < n; ++i) {
        double v = (Spec(i) - min_val) / norm_val; // 归一化到 [0, 1]

        uchar b = cv::saturate_cast<uchar>(255 * (1 - v));  // 蓝色通道
        uchar g = cv::saturate_cast<uchar>(255 * std::max(0.0, std::min(1.0, v * 2 - 0.5)));  // 绿色通道（中间增强）
        uchar r = cv::saturate_cast<uchar>(255 * v);        // 红色通道

        color_map.at<cv::Vec3b>(i / img.cols, i % img.cols) = cv::Vec3b(b, g, r);
    }
    return color_map;
}

Eigen::VectorXd ImageProcessor::denoiseSpecVector(const Eigen::VectorXd& Spec, const cv::Mat& img, int kernelSize) {
    int rows = img.rows;
    int cols = img.cols;
    cv::Mat SpecMat(rows, cols, CV_64F);
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            SpecMat.at<double>(i, j) = Spec(i * cols + j);
        }
    }

    cv::Mat denoisedMat(rows, cols, CV_64F);
    int halfKernel = kernelSize / 2;

    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            std::vector<double> neighborhood;
            for (int dy = -halfKernel; dy <= halfKernel; ++dy) {
                for (int dx = -halfKernel; dx <= halfKernel; ++dx) {
                    int ni = i + dy;
                    int nj = j + dx;
                    if (ni >= 0 && ni < rows && nj >= 0 && nj < cols) {
                        neighborhood.push_back(SpecMat.at<double>(ni, nj));
                    }
                }
            }
            std::sort(neighborhood.begin(), neighborhood.end());
            denoisedMat.at<double>(i, j) = neighborhood[neighborhood.size() / 2];
        }
    }

    Eigen::VectorXd denoised(rows * cols);
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            denoised(i * cols + j) = denoisedMat.at<double>(i, j);
        }
    }

    return denoised;
}

